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We study theoretically the far-from-equilibrium relaxation dynamics of spin spiral states in the 
three dimensional isotropic Heisenberg model. The investigated problem serves as an archetype for 
understanding quantum dynamics of isolated many-body systems in the vicinity of a spontaneously 
broken continuous symmetry. We present a field-theoretical formalism that systematically improves 
on mean-field for describing the real-time quantum dynamics of generic spin-1/2 systems. This is 
achieved by mapping spins to Majorana fermions followed by a 1/N expansion of the resulting two- 
particle irreducible (2PI) effective action. Our analysis reveals rich fluctuation-induced relaxation 
dynamics in the unitary evolution of spin spiral states. In particular, we find the sudden appearance 
of long-lived prethermalized plateaus with diverging lifetimes as the spiral winding is tuned toward 
the thermodynamically stable ferro- or antiferromagnetic phases. The emerging prethermalized 
states are characterized by different bosonic modes being thermally populated at different effective 
temperatures, and by a hierarchical relaxation process reminiscent of glassy systems. Spin-spin 
correlators found by solving the non-equilibrium Bethe-Salpeter equation provide further insight 
into the dynamic formation of correlations, the fate of unstable collective modes, and the emergence 
of fluctuation-dissipation relations. Our predictions can be verified experimentally using recent 
realizations of spin spiral states with ultracold atoms in a quantum gas microscope [S. Hild, et al. 
Phys. Rev. Lett. 113 , 147205 (2014)]. 

PACS numbers: 75.10.Jm, 05.40.-a 05.70.Ln, 


I. INTRODUCTION 

The equilibration of isolated quantum many-body sys¬ 
tems is a fundamental and ubiquitous question in physics. 
It plays a central role in understanding a broad range 
of phenomena, including the dynamics of the early uni¬ 
verse [1] , the evolution of neutron stars [2], pump-probe 
experiments in condensed matter systems [3] , and the op¬ 
eration of semiconductor devices [4]. The simplest per¬ 
spective on the problem is to recognize a dichotomy be¬ 
tween ergodic and non-ergodic systems. The former ex¬ 
hibit fast relaxation to local equilibrium states occurring 
at microscopic timescales, followed by a slower relaxation 
process to global thermal equilibrium described by classi¬ 
cal hydrodynamics of a few conserved quantities [5-7]. In 
contrast, non-ergodic systems possess an extensive set of 
conservation laws that prevent their thermalization [8, 9]. 

Recent theoretical and experimental investigations of 
strongly-correlated systems, however, suggest significant 
refinements to this dichotomy. For instance, certain 
systems can be trapped for long times in quasi-stationary 
“prethermalized” states with properties strikingly differ¬ 
ent from true thermal equilibrium [10]. Examples include 
nearly-integrable one-dimensional systems [11-15], and 
systems with vastly different microscopic energy scales 
in which slow dynamics results from the slow modes 
providing configurational disorder and thereby localizing 
the fast modes [16-18]. Even subtler examples of slow 


dynamics include the Griffiths phase of interacting 
disordered systems [19, 20] and translationally invariant 
systems in higher dimensions with emergent slow degrees 
of freedom [10, 21-29]. 

In this work, we discuss the emergence of slow dynam¬ 
ics and prethermalization in translationally invariant spin 
systems that possess continuous symmetries. In higher 
dimensions, these systems can exhibit thermodynami¬ 
cally stable symmetry broken phases along with gapless 
Goldstone modes. Here, we show that the relaxation dy¬ 
namics of low energy initial states that allow symme¬ 
try breaking upon thermalization is remarkably different 
from the relaxation of high energy states, thereby estab¬ 
lishing a connection between aspects of equilibrium and 
non-equilibrium phenomena in these system. In partic¬ 
ular, the slow Goldstone modes in the former case re¬ 
sult in the emergence of long-lived non-thermal states 
with a hierarchical relaxation dynamics that closely re¬ 
sembles aging in systems with quenched disorder. A 
non-perturbative treatment and beyond mean-field cor¬ 
rections are both found to be crucial for describing the 
relaxation process. 

We specifically study the dynamics of the three dimen¬ 
sional (3D) isotropic Heisenberg model initially prepared 
in a spiral state, see Fig. 1 (a). Our choice of spin spi¬ 
ral states is motivated by the following considerations. 
First, the winding of a spiral, Q, serves as a tuning pa¬ 
rameter for the energy density of the state. The full 
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FIG. 1. Relaxation of spin spiral states in the 3D 
isotropic Heisenberg model. (a) The system is pre¬ 
pared in a spin spiral state in the xy plane with the winding 
Q = (Q, Q, Q) as tuning parameter. The hgure illustrates the 
case Q = 7r/2. (b) The real-time evolution of the transverse 
magnetization M± for three different Q as indicated in the 
plot. For Q = Ttt/S, a hierarchical relaxation process emerges 
with a non-thermal plateau at intermediate times. The time 
scale is switched to logarithmic at tJ = 5 for better visibil¬ 
ity. (c) A global view of the spiral dynamics. Non-thermal 
plateaus appear near Q ~ 0, tt. 


spectrum of the Heisenberg model is traversed from ferro¬ 
magnetic (FM) to Neel antiferromagnetic (AFM) states 
upon sweeping Q from 0 to tt, respectively. In light of the 
eigenstate thermalization hypothesis (ETH) [30-32], the 
energy density of the initial state fully determines the fate 
of all local observables at late times in ergodic systems. 
One of the main objectives of this work is to understand 
the route toward thermalization of these states. Second, 
spiral states represent different mean-field solutions of 
the classical Heisenberg model, all of which are thermo¬ 
dynamically unstable with the exception of Q = 0,7r. 
The fluctuation-induced destruction of the initial order 
and the emergence of thermodynamically stable ordered 
or disordered phases at longer times is another question 
we address here. Lastly, spin spiral states have been re¬ 
cently realized in one and two dimensions using ultracold 
atoms in a quantum gas microscope [33, 34]. An exten¬ 
sion of these experiments to three dimensions makes a 
direct experimental scrutiny of our predictions possible. 

Our results indicate that spiral states tuned toward 
FM or AFM states exhibit a slow hierarchical relaxation 
and can come arbitrary close to a dynamical arrest, see 
Fig. 1 (b-c). Surprisingly, the relaxation dynamics is nei¬ 
ther compatible with the trivial relaxation to local ther¬ 


mal equilibrium and slow hydrodynamic evolution, since 
the spin current is not conserved, nor with the linearized 
dynamics of the collective modes, which predicts expo¬ 
nentially growing out-of-plane instabilities. In fact, we 
find the instabilities to self-regulate and slow down sig¬ 
nificantly. As we elaborate in the following sections, the 
physical phenomena discussed here are expected to gen¬ 
eralize to a broad range of models that exhibit a finite 
temperature phase transition between a disordered and 
a symmetry broken phase. 

The relaxation of the Neel spin spiral state with 
Q = TV has been previously studied in the ID Heisenberg 
model [33, 35-37]. In contrast to the 3D case studied 
in the present work, the ID Heisenberg model does not 
exhibit a symmetry broken thermal phase and in turn, 
cannot exhibit the type of prethermalization we discuss 
here. More recently, the dynamics of the Neel state in 
the Fermi-Hubbard model on an infinite dimensional 
Bethe lattice has been investigated [38], however, the 
approach to the steady state could not be studied due 
to the small effective exchange interaction. 

From a technical perspective, our investigation of the 
non-equilibrium dynamics of spiral states has been en¬ 
abled by developing a non-perturbative field theoretic 
formalism applicable to generic spin-1/2 systems for ar¬ 
bitrary initial states and geometries, which we refer to 
as the “Spin-2PI” formalism. This is achieved using 
a Majorana fermion representation of spin-1/2 opera¬ 
tors [39, 40], enlargement of the spin coordination num¬ 
ber by a replica-symmetric extension, and ultimately a 
systematic 1/N fluctuation expansion of the real-time 
two-particle irreducible (2PI) effective action [41, 42]. 

The recent rapid progress in the phenomenol¬ 
ogy of far-from-equilibrium quantum dynamics and 
its broad applications has been enabled by similar 
non-perturbative functional techniques. Examples 
include extensive studies of the 0{N) model in non¬ 
equilibrium [43-45], thermalization, prethermalization 
and non-thermal fixed points [10, 22, 46-48], particle 
production, reheating and defect generation in inflation¬ 
ary universe models [49-53], and dynamics of ultracold 
fermionic and bosonic gases [54-56]. The present work 
is the first to utilize this powerful technique to study the 
far-from-equilibrium dynamics of interacting quantum 
spin systems. 

Understanding the emergence of slow dynamics near 
thermodynamic phase transitions has implications 
reaching far beyond the domain of condensed matter 
physics. Eor instance, studies of non-equilibrium quan¬ 
tum fields in the context of inflation and early universe 
dynamics have suggested that the slowing down of 
quantum evolution near phase transitions is a plausible 
explanation for the large number of light particles and 
broken symmetries in the observable universe [57]. 
Given that the experimental verification of theories 
about early universe phenomena are typically rather 
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FIG. 2. The Spin-2PI formalism illustrated. The spins 
(green arrows) precess about a fluctuating exchange field (un¬ 
certain blue arrows). The quantum fluctuations of the ex¬ 
change field are mediated by a real vector bosons (wiggly 
lines) and are suppressed by a factor of 1/N, permitting a 
systematic expansion. 


indirect, experiments with synthetic many-body systems 
that allow precise monitoring of real-time dynamics 
close to phase transitions could play an important 
role in elucidating the emergence of slow evolution. 
The dynamics of various interacting spin systems have 
been already investigated in experiments with synthetic 
quantum matter, including domain formation in spinor 
condensates [58, 59], the precise measurement of the 
evolution of spin flips in the ground state of ID lattice 
spin systems [60-63], quantum coherences in long range 
models [64], and the relaxation dynamics of spin spiral 
states in ID and 2D Heisenberg models [33, 34]. The 
experimental observation of the dynamical phenomena 
discussed here are thus expected to be within close reach. 

This paper is organized as follows: In Sec. II, we intro¬ 
duce the Spin-2PI formalism, a technique we develop to 
study the dynamics of interacting spin systems. Comple¬ 
mentary technical details are presented in App. A. We 
discuss the relaxation of spin spiral states in Sec. III. 
The phenomenon of dynamical slowing down and arrest 
will be presented Sec. Ill A, the long-time thermalization 
in Sec. IIIB, and the dynamic formation of correlations 
and instabilities in Sec. IIIC. We conclude our findings 
in Sec. IV. 


II. THE SPIN-2PI FORMALISM 

Consider a generic Hamiltonian describing the pairwise 
interaction between localized spin degrees of freedom on 


a given lattice L: 

^ = ( 1 ) 

j,ke\- 

where V is an arbitrary interaction, j and k denote lat¬ 
tice sites, and are spin-1/2 operators. Summation 

over the repeated spin indices is assumed. Had S been 
classical angular momentum variables, the Hamiltonian 
dynamics of the system would be governed by the (non¬ 
linear) Bloch equation: 

( 2 ) 

kei 

In case of quantum spins, the Bloch equation only de¬ 
scribes the evolution of the spin expectation values (S) 
to the extent of which the mean-field Ansatz {SjSk) ~ 
{Sj) {Sk) is valid. The latter, however, is only justified for 
lattices with large coordination number, high spin parti¬ 
cles, or in the presence of a high-temperature bath. The 
crucial role of quantum fluctuations in the dynamics of 
isolated spin-1/2 systems in finite dimensional lattices is 
beyond the reach of semi-classical methods, and demands 
a more careful treatment. 

Here, we propose a formalism for transcending the 
mean-field approximation for spin evolution by a system¬ 
atic inclusion of quantum corrections. This is achieved 
using functional methods and a variant of the large-V 
expansion technique. As a first step, we construct an 
auxiliary model in which each spin is replicated N times, 
and each bond is promoted to bonds between the 
replicas, with equal weight but with an overall scale fac¬ 
tor of 1/N. The Hamiltonian of the auxiliary model is 
written as: 

*»=5 s (s t- <2) 

jEL \/cGL cr' = l / cr = l 

The initial state |lFo) is also subsequently promoted to an 
uncorrelated product in the replica space, 0^=1 l^o)cr- 
The original problem is recovered by setting N = 1. We 
refer to the sum appearing in the parentheses in Eq. (3) as 
the exchange field operator, (pj^ which plays the role of an 
effective fluctuating magnetic field with which the spins 
interact. The described large-V construction effectively 
increases the coordination number of each spin, 2 :, to Nz^ 
thereby suppressing the fluctuations of (pj according to 

the law of large numbers, p^ = p^ j + 0{1/\/Nz)^ where 
p^ - = {p-) is the mean exchange field. In the limit of in¬ 
finite V, the exchange field operator becomes effectively 
classical such that mean-field dynamics of the original 
model H emerges as the asymptotically exact description 
of the dynamics in limiv^oo For large but finite V, 
the fluctuations of p are small but not negligible, and can 
be systematically incorporated into the dynamics order 
by order in 1/N. This program can be carried out within 
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the functional method of two-particle irreducible (2PI) 
effective actions. Crucially, truncating the expansion at 
a finite order in 1/N and taking the limit N ^ 1 yields 
non-perturbative and conserving approximations for the 
spin dynamics. We refer to this method as the Spin-2PI 
formalism, which is illustrated schematically in Fig. 2. In 
brief, spins precess about a self-consistently determined 
exchange mean field, and quantum spin fluctuations are 
mediated by the local and non-local exchange of a real 
vector boson whose propagator is suppressed by a factor 
of 1/A^. 

In the remainder of this section, we briefly outline the 
field theoretical developments that underlie the Spin-2PI 
formalism. Complementary technical details are given 
in App. A. A path integral for the spin-1/2 operators 
is constructed using a representation invoking Major ana 
fermions [39, 40]: 


to Eq. (4), the local spin expectation value is proportional 
to the fermion tadpole: 

= (7) 

We obtain the real-time evolution equations for 0, P and 
using the 2PI effective action formalism [41]. The 
effective action r[0, P, cp^] is found by sourcing Q, V and 
and performing Legendre transformations: 

- \ tr[Do ^r>] + (8a) 

g] + ^(pj)Q^ip^ + r2[G,v], 

(8b) 


Si =X 77^.. (4) 

The Majorana operators at each site {tj] satisfy 

the Clifford algebra = SjkS^^, from which the 

SU{2) algebra for spins [S^,S^] = iSjk £ai 3 'y Sj and the 
Casimir condition S| = 3/4 follow. The latter ensures a 
faithful spin-1/2 representation without introducing any 
unphysical states and obviates the necessity of using con¬ 
straint gauge fields in contrast to the Schwinger slave 
particle approach [65]; see Appendices of Ref. [66] for 
a detailed treatment of the Majorana representation for 
spin-1/2 operators. 

Replacing the spin operators in H using Eq. (4), the 
Hamiltonian is mapped to that of a many-body system of 
Majorana fermions with quartic interactions. The large- 
N program can be identically followed by replicating the 
slave Majorana particles and assigning a replica index 
to each. We proceed by constructing a path integral for 
the Majorana fermions using fermionic coherent states on 
the closed time path (CTP) Schwinger-Keldysh contour. 
The Lagrangian is given as: 

iGL cr=l i,/cGL cTi,0-2 = 1 

iVj X VjT’''" iVk X Vkf (5) 

The exchange field is introduced by a Hubbard- 
Stratonovich decoupling of the quartic term using a real 
vector boson cpj on each lattice site. The non-equilibrium 
exchange mean field exchange field fluctuation prop¬ 
agator P, and the Majorana propagator Q are introduced 
as: 


<^c(l) = (<^(1)), 

iV{l, 2) = {Tc[<pil) <pi2)]) - (^,(1) (^,(2), 

ig{l,2) = {Tc[v{l)vm)- (6) 

The integer variables are shorthand for the bundle of lat¬ 
tice site, contour time, spin and replica index. According 


The bare Majorana and exchange propagators are 
given as = idt^S{l,2) and 'Dq^{1,2) = 

N (V~^) respectively. M[cp^] is the leading or¬ 

der (LO) self-energy (see Eq. A9). The evolution equa¬ 
tions follow from making f stationary with respect to G, 
P, and (p^, see Eqs. (A7a)-(A7c). 

Save for r 2 [G^'^]^ the rest of the terms appearing 
in r[(^^,0,P] scale as 0{N) and together comprise the 
leading-order (LO) approximation. The next-to-leading- 
order (NLO) corrections and beyond are represented by 
r2[0,T>] which formally corresponds to the sum of 2PI 
vacuum diagrams arising from the cubic interaction ver¬ 
tex: 


>^int [17 : ^] 


% % 

-(p- iV'XVV = ^ 



nl3',o- 


(9) 


The 1/N expansion of Tint to the next-to-next-to-leading 
(NNLO) order is diagrammatically given as: 


^int[G,V] = Q--Q 

LO ~ OiN) 



NLO ~ e>(l) 



NNLO ~ e>(l/Ar) 


( 10 ) 

We have used the stationarity condition, Eq. (A7c), to 
omit (p^ in favor of G in the LO interaction terms. The 
Eeynman diagram rules are given in Sec. A 2. 

Truncating the 1/N expansion of f at a finite order 
and setting N = 1 yields systematic improvements of 
the mean-field spin dynamics. The ensuing approximate 
theories are self-consistent and non-perturbative by con¬ 
struction, and respect the conservation laws associated 
to the global symmetries of the microscopic action, such 
as magnetization and energy. The latter is crucial for the 
long-time stability of the non-equilibrium dynamics. 

The Bloch equation is recovered upon truncating f at 
the LO level, see Sec. A3. Truncations at NLO and be¬ 
yond give rise to memory effects due to the dynamical 
fluctuations of the exchange field and result in a two- 
time Kadanoff-Baym integro-differential equation instead 
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of the mean field Bloch equation, see Eqs. (Al3a)-(A14b). 
Finally, higher order correlators, in particular the spin- 
spin correlator ix(l,2) = {Tc[S(1)S{2)]) — {S{1)) (^(2)), 
can be reconstructed with the knowledge of Q and V by 
solving the non-equilibrium Bethe-Salpeter integral equa¬ 
tion on the Schwinger-Keldysh contour, see App. A3. 

We remark that in systems with large spin coordina¬ 
tion number 2 ;, fluctuations of the exchange held are in¬ 
herently suppressed and the expansion parameter is more 
accurately identified with l/{zN). Therefore, the large- 
N expansion of the Spin-2PI effective action in models 
with ^ ^ 1 is expected to be controlled and rapidly con¬ 
verging, even after taking the limit N ^ 1. Studies of 
the 0{N) model show that the most important correc¬ 
tion to the mean-field (LO) approximation is captured 
by the NLO “fluctuation-exchange” diagram, along with 
negligible quantitative corrections from the subleading 
terms [67, 68]. 

The replica-based 1/N expansion proposed here differs 
from the usual semi-classical 1/S expansion in significant 
ways even though they improve upon the same mean- 
field limit. For instance, the replicated Fock space of a 
single spin is reducible and has many more states com¬ 
pared to a pure spin-A^/2 representation. A technical 
advantage of our approach is that it preserves the under¬ 
lying spin-1/2 degrees of freedom, which in conjunction 
to the Majorana representation leads to the familiar dia¬ 
grammatic and functional methods. As discussed before, 
these tools significantly simplify and streamline the cal¬ 
culation of higher order corrections. Furthermore, there 
is no preferred axis for spin quantization in the Spin-2PI 
formalism, allowing us to study magnetically ordered and 
disordered states in a unified way. 


III. RELAXATION OF SPIN SPIRAL STATES 
IN THE 3D HEISENBERG MODEL 

In this section, we investigate the unitary evolution of 
the spin spiral state on a 3D cubic lattice, 

|sp(Q))=e-*^iQ'^^'^I (g)K)^., (11) 

under the isotropic Heisenberg Hamiltonian H = 
• Sj using the Spin-2PI formalism devel¬ 
oped in the previous section. Here, \^)j denotes the 
x-polarized state on lattice site j. The spiral is prepared 
in the x^-plane with a winding wavevector Q. We as¬ 
sume ferromagnetic couplings J > 0 for concreteness, 
even though the sign of the J does not affect the uni¬ 
tary evolution due to the time reversal symmetry of the 
Heisenberg model. 

The spiral state |sp(Q)) is a simultaneous eigenstate 

of 5a(Q) = fa f-ziQa), a = x,y, z, where fa and TZziQa) 
denote the translation by one lattice site along the a-axis 
and rotation by angle Qa about the 2 ;-axis, respectively. 
The translation and rotation symmetries of the isotropic 


Heisenberg model imply [H^Sa] = 0, so that the spi¬ 
ral state |sp(Q)) remains a simultaneous eigenstate of 
5a (Q) at all times in the course of unitary evolution. As 
a result, the out-of-plane magnetization (5|(t)) vanishes 
identically, and the spiral magnetic order with the initial 
winding Q persists at all times. The transverse magne¬ 
tization, 

jei 

is the only degree of freedom at the level of single spin 
observables. Also, Mx(k,t) = 0 for k 7^ Q. We remark 
that even though the magnetization dynamics is signifi¬ 
cantly constrained at the level of single spin observables 
by symmetries, arbitrary spin correlations are allowed to 
form in the course of evolution, including both in- and 
out-of-plane spin correlations at arbitrary wavevectors. 

A simplifying aspect of the present problem is that 
the apparently broken translation symmetry of the spi¬ 
ral state can be restored using an “unwinding” unitary 
transformation Uq = under which the spi¬ 

ral state transforms into a uniform x-polarized prod¬ 
uct state |!^o) = Uq |sp(Q)) = 0^- |^)j- The unwind¬ 
ing transformation, however, transforms the Hamiltonian 
H ^ H = UqHUq to an anisotropic Heisenberg model 
with a Dzyaloshinskii-Moriya term: 


H = -JJ2[S]SI+ cosQ• (R,- -Rfe) + Spl) 




(13) 


The translation invariance of the initial state in the 
spiral frame significantly simplifies the structure of the 
Spin-2PI equations: Q and E become local in the real 
space while V depends only on the distance between the 
sites. These simplifications hold for arbitrary truncations 
of Tint- Additionally, the bosonic self-energy 77 becomes 
local in the real space at the NLO truncation. The mag¬ 
netization is non-vanishing only along the x-direction 
in the spiral frame due to the symmetry considerations 
mentioned earlier. The quantities calculated in the 
spiral frame can be readily transformed to the lab frame 
using appropriate rotations. In particular, Eq. (7) gives 
M^(Q,t) = (1/2) t) with Q calculated in the 

spiral frame. We choose the winding to be along the 
diagonal direction Q = (Q, Q, Q) hereafter and refer to 
the spiral winding with the single scalar Q G [0,7r]. 


At the LO level, the spin dynamics is governed by the 
Bloch equation, Eq. (2). The exchange mean held cp^ 
is parallel to the local magnetization at all lattice sites 
in a spiral state, implying the absence of any dynamics. 
In other words, the spiral states are fixed points of the 
mean-field dynamics for all windings Q. 

Going beyond the LO dynamics and including the 
exchange held fluctuations by taking into account the 
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NLO corrections, the spiral state exhibits an intrigu¬ 
ing fluctuation-induced relaxation dynamics. States with 
different windings have different energy densities, along 
with different strength of in-plane and out-of-plane spin 
fluctuations, and are found to relax in strikingly differ¬ 
ent ways. As we discuss below, these factors conspire to 
give rise to a non-trivial hierarchical relaxation scenario 
for spiral states lying close to thermodynamically stable 
orders, exhibiting prethermalization [10], and dynamical 
arrest resembling glassy systems [69]. 


A. Relaxation and dynamical arrest of the 
transverse magnetization 

The spiral state for Q = 0 is a fully polarized FM 
eigenstate of the Heisenberg model and is therefore 
stationary. The Q = tt spiral, on the other hand, 
corresponds to an uncorrelated Neel state which in three 
dimensions has a large overlap with the correlated AFM 
state lying at the upper end of the spectrum of the FM 
Heisenberg model. As a result, the system is expected 
to achieve a steady state marked with a finite staggered 
magnetization after a short course of dephasing dynam¬ 
ics, provided that the generated effective temperature 
is below the ordering temperature. The evolution of 
is shown in Fig. 1 (b) for several choices of Q, along 
with a global surface plot for Q G [0, tt] and tJ G [0, 30] 
in Fig. 1 (c). The stationarity of the FM state {Q = 0) 
and the rapid settlement of Neel state {Q = tt) to 
a steady state with finite staggered magnetization is 
observed. 

Short-time dephasing dynamies — For all Q, the first 
stage of dynamics is a short-time relaxation of the form 
M± ^ 1/2 — arising from the dephasing between 
the eigenstates that overlap with the spiral. A straight¬ 
forward calculation using the short-time expansion 

im) = {S)o + it {[H,S])o + +■■■ 

gives ^ (cosQ — 1)^- The values of z/g extracted 

from the numerically obtained M±_ are in agreement 
with the exact result, see Fig. 6. The second stage of 
relaxation dynamics depends on the winding of spiral 
and is either directly thermalizing, or exhibits long-lived 
prethermalized states preceding the true thermalization. 
We discuss these cases separately. 

Spiral states with Q ^ 7r/2— Spin spiral states with 
Q ^ 7T12 have a high energy density with respect to 
both the FM and the AFM state. Thus, Q ^ 'k j2 spiral 
states overlap with a large number of eigenstates of the 
Heisenberg model. Such a broad superposition of states 
lead to fast dephasing which is found to be within a 
few exchange times. Our results indicate a rapid onset 
of exponential decay ^ e~^^^ with the fastest rate 
occurring at Q = 0.55(l)7r ~ 7r/2. 

Spiral states with Q ^ 0 and Q ^ ir — A complex multi- 
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FIG. 3. Thermalization of the spin spiral state, (a) The 

effective inverse temperature of local spin Tgpin and local ex¬ 
change field fluctuations Tfluct. obtained from the fluctuation- 
dissipation relations in steady state. The two temperatures 
are in agreement for spiral windings near Q ~ 7r/2, support¬ 
ing the true thermalization of the system. The temperatures 
calculated in the prethermalized plateaus Q ~ 0, tt (shaded re¬ 
gions) disagree with each other, and generically differ from the 
temperature of the true thermal states that emerge at later 
times. Inset: the temperature ksT as a function of Q (same 
data as in the main panel) displays a resonance from positive 
infinite temperature to negative infinite temperature at the 
classical duality point Q = 7r/2. (b) The approach of Tfluct. 
to steady state (light to dark) as obtained from fluctuation- 
dissipation relations for Q = 7r/4 (left) and Q — ti (right). 
The steady state temperatures are shown on the plots. 


scale relaxation scenario emerges for spirals with wind¬ 
ings tuned to Q ^ 0 and Q ^ tt, lying close to FM and 
AFM magnetic orders, respectively. The transverse mag¬ 
netization exhibits an intermediate plateau for these ini¬ 
tial states which appears continuously upon tuning Q, see 
Fig. 1. The plot of M±_ shown in Fig. 1 (b) for Q = Ttt/S 
displays the intermediate plateau followed by relaxation 
at later times. As Q is tuned closer toward 0 or tt, the 
lifetime of the plateau increases abruptly and the mag¬ 
netization comes to a dynamical arrest. We investigate 
the nature of such long-lived plateaus in more detail in 
the following sections. 
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FIG. 4. The evolution of spin correlations. Top panels: Growth rate of out-of-plane instable modes obtained from a linear 
response analysis. Bottom panels: Numerically calculated correlation function (*§k*§-k)(^) = as a function of the 

lattice wavevector k = (/c, /c, k) within the Spin-2PI formalism including NLO corrections, (a) Q — Stt/S, (b) Q — 7r/2, and (c) 
Q = 37r/4. The inset in (c) shows the connected part of the in-plane correlations *§Zk)(^)- 


B. Prethermalization vs. Thermalization 

Due to the non-integrability of the 3D Heisenberg 
model, the energy distribution of spin fluctuations is ex¬ 
pected to approach a thermal population in the long time 
limit, according to the eigenstate thermalization hypoth¬ 
esis (ETH) [30-32]. We investigate the nature of steady 
states emerging in the dynamics by calculating the spin- 
spin correlation and response functions, corresponding 
to the Keldysh (K) and retarded (R) components of the 
CTP spin-spin correlator x(t,t'), by solving the non¬ 
equilibrium Bethe-Salpeter equation (see App. A3). At 
thermal equilibrium, these quantities are related via the 
bosonic fluctuation-dissipation relation (FDR): 

iX^{^) = —2 coth(u;/2/cBT) Im[x^(ci;)] , (14) 

where T is the effective temperature. Here, uj refers to the 
Fourier frequency in the time difference t—t' in the steady 
state achieved at long times. Likewise, one can define an 
effective temperature for the exchange field fluctuations 
using the bosonic FDR between and V^. We re¬ 
fer to the temperatures obtained from local spin % and 
exchange fluctuations V as Tgpin and T^uct.i respectively. 

The effective temperatures obtained from the FDR in 
the steady state are shown in Fig. 3(a). For all spiral 
windings Q, we find that FDR is satisfied to an excellent 
degree for both local spin and exchange fluctuation cor¬ 
relators once the steady state is reached, see Fig. 3(b). 
However, as we discuss below, the effective temperature 
obtained from spin and exchange fluctuations may 
disagree with each other. This allows us to distinguish 
prethermalization from true thermalization. 

Thermalization of spiral states with Q ^ 7r/2— For a 
range of spiral wavevectors 7r/4<(5^37r/4, the steady 
state temperatures obtained from all bosonic modes, i.e. 
local and non-local in- and out-of-plane spin and ex¬ 
change field fluctuations, agree with each other, suggest¬ 


ing the complete thermalization of the system and in ac¬ 
cordance with the ETH. 

Spirals with Q = nj2 flow to an infinite temperature 
thermal state, which is understood from the duality 
Q ^ TT — Q, J ^ —J present in the classical Heisenberg 
model. This classical duality extends to the quantum 
Heisenberg model in the high temperature regime. The 
duality point Q = 7t/2 further marks the resonance 
from positive temperatures for Q < 7r/2 to negative 
temperatures for Q > 7r/2, see the inset of Fig. 3 (a). The 
T < 0 thermal states of the FM Heisenberg model with 
coupling —\J\ corresponds to T > 0 states of the AFM 
Heisenberg model with coupling +|J|, and vice versa. 
Negative temperature states naturally arise in isolated 
systems with bounded energy spectra as legitimate 
thermal states and occur when the initial energy density 
lies closer to the upper edge of the energy spectrum. 

Prethermalization of spiral states with Q ^ 0,7r— For 
spiral states with Q ^ 0, tt, where the system develops a 
prethermal plateau, the effective spin and exchange field 
fluctuation temperatures disagree, even though the FDR 
is satisfied well for each mode individually. This find¬ 
ing supports the prethermalized nature of such steady 
states. It is understood that the temperatures calcu¬ 
lated within the prethermal plateau [shown as shaded 
regions in Fig. 3 (b)] correspond to the effective temper¬ 
ature of individual modes, and not the true thermody¬ 
namical temperature. We expect the two temperatures 
to approach each other at longer times once the system 
exits the prethermalized plateau and progresses toward 
a fully thermalized state. 

The spiral state with Q = 0 is an exact ground state of 
the system and FDR yields T = 0 as expected. In con¬ 
trast, the Q = 7T state approaches a finite temperature, 
which is understood by the fact that the uncorrelated 
Neel state must be “dressed” with spin correlations 
before the steady state is reached, see inset of Fig. 3 (a). 





















The disparity between the evolution of Q = 0 and Q = tt 
states reveals the quantum mechanical nature of spins, 
and the breakdown of the classical duality Q ^ ir — Q 
in the low temperature regime. 

The 3D Heisenberg model exhibits a finite temperature 
equilibrium phase transition from the disordered para¬ 
magnetic phase to the ordered FM or AFM phase, de¬ 
pending on the sign of the exchange coupling J. For the 
spiral at Q = TT, the FDR of the spin fluctuations are 
not well fulfilled at accessible times while those for ex¬ 
change field fluctuations are. The temperature extracted 
from the latter |Tfluct.(Q = '7r)| = 0.82J lies below the 
the AFM ordering temperature = 0.946(1) J. The 

latter has been obtained from quantum Monte Carlo sim¬ 
ulations [70]. 

Crucially, the near-thermal distribution of fluctuations 
in the prethermalized plateaus, the stability of FM/AFM 
ordered phases at finite energy densities in the 3D Heisen¬ 
berg model, and the proximity of Q ^ 0, tt spiral states to 
these stable orders allow us to draw a connection between 
the long-time stability of such spiral states and spon¬ 
taneous symmetry breaking at equilibrium: the spiral 
winding Q sets the energy density of the system and sub¬ 
sequently the effective temperature T{Q) in the prether- 
mal state. T{Q) approximately dictates the magnitude 
of spin fluctuations on the top of the spiral states which 
locally resemble either FM or AFM for Q ^ 0,7r. De¬ 
pending on Q, T{Q) can either he below or above the 
critical transition temperature, T™ or thereby 

providing an approximate condition for the local stabil¬ 
ity of the spiral order. We will study the global instability 
of the spiral states and their destruction at longer times 
in the next section. 

According to the above discussions, the connection 
made between emergence of slow dynamics and sym¬ 
metry breaking at equilibrium essentially hinges on the 
eigenstate thermalization hypothesis and the Mermin- 
Wagner theorem. Therefore, this connection is expected 
to reach beyond the present discussion, and to general¬ 
ize to a broader range of initial states and models that 
exhibit spontaneous continuous symmetry breaking. 


C. Instabilities and Correlations 

The dynamical stabilization of spirals near the FM 
and AFM orders, and consequently the appearance of 
prethermal plateaus, was understood on the basis of 
thermodynamical arguments in the previous section. 
However, even though the spiral states are fixed points 
of the mean-field dynamical equations, they are unstable 
and have a tendency to form out-of-plane textures as the 
energy of spiral states can be reduced by an appropriate 
out-of-plane tilt. Therefore, in a thermodynamical 
ensemble where arbitrary out-of-plane fluctuations are 
allowed, these saddle points fail to give rise to symmetry 
broken states, leaving Q = 0 and Q = tt as the only 


thermodynamically stable orders in the Heisenberg 
model. Therefore, the present situation must be re¬ 
garded from the perspective of quantum dynamics, i.e. 
the unitary evolution of a pure state |sp(Q)) rather than 
the statistical fluctuations in a mixed thermodynamical 
ensemble. Here, the system remains in a pure state 
at all times and the magnetic order is confined to the 
xy spiral plane due to the symmetries discussed at 
the beginning of Sec. HI. It is therefore conceivable 
that symmetry-protected dynamical constraints allow 
thermodynamically unstable saddle points to become 
long-lived states in the course of unitary dynamics. 

Out-of-plane instability .— The out-of-plane instability of 
the spiral state in the Heisenberg model can be studied 
either by performing a linear response analysis of the 
Bloch equations, or similarly from the Holstein-Primakoff 
spin-wave analysis. Either way, the dispersion of out-of- 
plane spin-waves forming on the top of the spiral is found 
as cjk = where: 

3 

ek = —JS ^ [(1 + cos Q • Bd) cos k • — 2 cos Q • , 

d=l 

3 

Z\k = — JS ^ [(1 — cos Q • Bd) cos k • . (15) 

d=l 

Unstable modes arise when cjk assumes imaginary val¬ 
ues. Except for Q = 0,7r/2,7r, one always finds such 
unstable modes: for Q < 1^12^ the fastest growing mode 
is k = {k,k,k) with k = cos“^[cos^((5/2)] along with a 
sharp cutoff \k\ < Q; for Q > 7r/2, unstable modes occur 
for \k — 7r\ < Q, with the fastest mode always being the 
staggered k = ir mode, independently of Q. 

A simple estimate for the lifetime of the prethermal 
plateaus is obtained by calculating the time it takes for 
typical unstable out-of-plane collective mode to grow 
to 0(1). The rationale behind this estimate is that 
the in-plane order can not coexist with strong enough 
out-of-plane fluctuations. Expanding around Q = 0,7r, 
we obtain a scaling t ^ 1/Q^ for the EM-like and 
t ^ l/(7r — Q) for the AEM-like spirals, up to logarithmic 
corrections. However, the lifetime of plateaus as found 
from the Spin-2PI formalism exceeds the above esti¬ 
mates; in particular, as Q is tuned closer toward 0 or tt, 
we observe a faster increase of the lifetime of prethermal 
plateaus. As we discuss below, the increased lifetime 
can be explained on the basis of the self-regulation of 
out-of-plane spin fluctuations. 

The top panels in Eig. 4 show Im[co’k] for several values 
of Q, along with the evolution of equal-time out-of-plane 
spin correlations ^Xk^’^(t,t) = {S^{t) Sf^{t)) obtained 
by solving the non-equilibrium Bethe-Salpeter equation 
in the NLO approximation, bottom panels. Eurther, the 
connected part of in-plane correlations = 

{S^(t) SZ\^{t)) — (^k (^))(^-k(^)) shown in the inset 
of panel (c). 
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FIG. 5. Dynamics of exchange field fiuctuations. (a) The out-of-plane (top) and in-plane (bottom) exchange field 
fluctuations as a function of time and momentum k = (/c, /c, k) for Q = 7r/4 (left) and Q = Ttt/S (right). The red lines indicate 
the most enhanced mode in the long time limit; the blue dashed line in the lower right plot corresponds to the k = Q in-plane 
mode which initially exhibits the strongest enhancement of correlations, (b) The evolution of the late-time most enhanced 
mode for Q = (left) and Q — (right). In cases where the system thermalizes, left column, SU{2) symmetry 

emerges in the long time limit, while it is broken in the prethermal case Q — Ttt/S and for Q = tt, right column. In the latter 
case, the system can exhibit true long-range order provided its effective temperature is below the critical temperature of the 
equilibrium phase transition and thus be thermal and simultaneously break SU{2) symmetry. 


At t = 0, spin correlations are zero in accordance 
with the initial spiral state |sp(Q)) being an uncorre¬ 
lated product state. The out-of-plane correlations form 
at times t ^ J. The most enhanced correlations coin¬ 
cide with the wavevector predicted by the linear response 
analysis to a good degree. The sharp cutoffs predicted by 
this analysis are found to be smeared, which is expected 
due to the mode coupling embedded in our self-consistent 
approach. The time scale for the formation of correla¬ 
tions is found to be on the order of the dephasing time, 
reflecting the fact that the short-time dephasing dynam¬ 
ics and formation of correlations are manifestations of 
the same phenomenon. 

For spiral states that thermalize within the numeri¬ 
cally achievable time scales, we observe a smooth shift in 
both in-plane and out-of-plane spin correlations from the 
initial Q-dependent enhanced modes to either k = 0 or 
/c = TT, depending on whether Q < 7r/2 or Q > 7r/2, 
respectively [see Fig. 4(a), (c), and the inset]. Even 
though the linear response analysis correctly indicates 
the wavevector of the fastest growing out-of-plane mode, 
the spin correlations rapidly saturate to their maximum 
values, as opposed to an unbounded exponential growth. 
A similar rapid dynamical regulation of the growth of un¬ 
stable modes was previously reported in Ref. [22, 49] in 
the context of parametric resonance in the 0{N) model. 
In the present context, this phenomenon explains why 
the lifetime of the plateaus exceeds the estimate obtained 
from the linear response analysis, and indicates the im¬ 
portant role of mode coupling between spin waves and 
the necessity of non-perturbative treatments. 

For spiral states that exhibit long-lived prethermal 
plateaus, we study the exchange field correlations P, a 


quantity that is closely related to x but can be calculated 
for much longer times with less computational resources. 
The evolution of t) and t) for Q = 7r/4 

and Q = Ttt/S are shown in Fig. 5 (a). The former corre¬ 
sponds to a spiral state that thermalizes at about 20 
while the latter exhibits a magnetization plateau up to 
tm ^ 20 as shown in Fig. 1. For t < tm, the most 
enhanced in-plane mode occurs ai k = Q which upon 
demagnetization smoothly switches to /c = tt for t ^ tm- 
The most unstable out-of-plane mode is always at k = tt. 

As a further check for thermalizing behavior, we study 
the restoration of the SU (2) symmetry in the exchange 
field fluctuations {t^t) (l/2)P^~’^(t, t), see 

Fig. 5 (b). For Q = Tr/4 and Q = Tr/2, we find that the 
SU (2) symmetry is restored at longer times (left column), 
while for Q = Ttt/S and the Neel initial state Q = tt, the 
SU{2) symmetry remains broken at all accessible times 
(right column). Notably, the out-of-plane fluctuations for 
Q = Ttt/S is found to be an order of magnitude stronger 
than the Q = tt, in agreement with the previously men¬ 
tioned existence of an unstable out-of-plane mode for the 
former state and its absence in the latter. 

The magnitude of in-plane fluctuations remain es¬ 
sentially constant in the plateau for Q = Ttt/S (top 
right) while the out-of-plane fluctuations monotonically 
increase and reach a maximum at t ^ 20 ^ tm, pre¬ 

cisely when the prethermal magnetization decays. This 
finding connects the decay of the the spiral to the growth 
of out-of-plane fluctuations. The time tm also marks a 
reversal in the trend of out-of-plane and in-plane corre¬ 
lations. Even though this change indicates a first step 
toward establishing 5'T/(2)-symmetric correlations, the 
condition is far from being satisfied at t ^ tm and is 
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FIG. 6. Comparison between Spin-2PI and semi- 
classical dynamics. The Spin-2PI results (solid black lines) 
are compared to the semi-classical dynamics obtained from 
the dTWA (dashed blue lines). The short time analytic result 
from Sec. Ill A is also shown for reference (thick red lines). 
The long time dynamics of the two methods are significantly 
different. In particular, the dTWA is not capable of describing 
the long-lived prethermal plateaus in contrast to the Spin-2PI 
formalism. The scale of the time axis is switched from linear 
to logarithmic at tJ = 5 for better visibility. 





bound to occur on much longer time scales, indicating 
a hierarchical relaxation scenario with the relaxation of 
magnetization preceding the relaxation of correlations. 

The appearance of long-lived prethermal states and 
the hierarchical relaxation is reminiscent of aging dy¬ 
namics in classical structural glass models with quenched 
disorder [69] and kinematically constrained models [72]. 
Similar multi-scale glassy relaxation dynamics has been 
recently reported in the quench dynamics of fermions in a 
nearly-integrable ID model using a different method [15]. 

Comparison to semi-classical methods — According to 
the discussions presented so far, the relaxation dynam¬ 
ics of the spiral state accompanies the formation of 
in- and out-of-plane quantum correlations in the sys¬ 
tem. In order to study the role of correlations fur¬ 
ther, we compare our predictions with the results ob¬ 
tained from the discrete truncated Wigner approxima¬ 
tion (dTWA) [73, 74], a variant of the semi-classical TWA 
method [75] that relies on mean-field trajectories. The 
magnetization obtained from Spin-2PI (solid black lines) 
and dTWA (dashed blue lines) are compared in Fig. 6. 
The two methods generically agree with the analytic 
short time expansion (red thick lines), with the excep¬ 
tion that dTWA does not reproduce the correct short 
time dynamics for small Q, cf. Q = 7r/4 in (a) [76]. The 
two methods, however, predict strikingly different long 
time dynamics. Even though dTWA exhibits some de¬ 
gree of dynamical slowing down for FM-like and AFM- 
like spirals, it produces neither the prethermal plateau 
for Q = 77r/8, nor the finite steady-state magnetization 
for Q = TT. We note that the latter is supported by exact 


QMC calculations. 


IV. CONCLUSIONS AND OUTLOOK 

We formulated a non-perturbative and conserving 
field theoretic technique for describing the far-from- 
equilibrium quantum dynamics of strongly interacting 
spin-1/2 systems for arbitrary lattices and initial states. 
Referred to as the Spin-2PI formalism, this method sys¬ 
tematically improves upon the mean-field description by 
including quantum fluctuations by means of an asymp¬ 
totic 1/N expansion, which is controlled in models with 
intrinsically large lattice coordination number. 

We utilized the Spin-2PI technique to study far-from- 
equilibrium phenomena in spin systems with continuous 
symmetries. Specifically, we explored the relaxation dy¬ 
namics of spin spiral states in the 3D Heisenberg model, 
treating the spiral winding Q as a tuning parameter. Go¬ 
ing beyond the trivial mean-field (LO) dynamics by in¬ 
cluding the NLO correction, we found the spiral states 
with different windings to relax in remarkably different 
ways. In particular, spiral states resembling FM and 
AFM ordered states, corresponding to Q ^ 0 and tt 
respectively, get trapped for long times in non-thermal 
states, i.e. “false vacuums” whose lifetime diverge as the 
windings are tuned to Q = 0 or tt. In contrast, the spiral 
states far from Q = 0, tt relax rapidly. 

We calculated the effective temperature of spin and ex¬ 
change field fluctuations from the fluctuation-dissipation 
relation. For Q ^ Trj2 spiral states, all modes reach a 
single temperature, supporting full thermalization in ac¬ 
cordance with the eigenstate thermalization hypothesis. 
In contrast, the different bosonic modes of prethermaliz- 
ing spirals settle at different temperatures. 

We investigated the dynamical formation of corre¬ 
lations and found that the collective modes predicted 
to be unstable from a linear response analysis, are 
self-regularized at rather short time scales, demon¬ 
strating the importance of the nonlinear effects and 
non-perturbative treatments. The growth of out-of-plane 
fluctuations cause the eventual decay of the prethermal 
states. The restoration of SU (2) symmetry occurs much 
later after the decay of magnetization, suggesting a 
hierarchical relaxation reminiscent of coarsening and 
aging in classical glassy systems. Our results can be 
tested readily in ultracold atoms experiments with two- 
component Mott insulators in 3D optical lattices, such 
as a 3D extension of the experiments in Refs. [33, 34]. 

This work can be extended in several directions. A 
straightforward extension is to investigate the relaxation 
of spiral states in anisotropic models, or in lower di¬ 
mensions. Another immediately accessible direction is to 
study spin systems with long-range interactions, as real¬ 
ized for instance with Rydberg atoms, polar molecules, or 
trapped ions, and their instability toward dynamic crys¬ 
tallization. The effect of small deviations from the initial 
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spiral state on the quantum evolution could be studied as 
well. We expect our predictions to carry over to the case 
of weakly disordered initial states, provided that the de¬ 
viations from the pristine spiral state remain small by the 
time the prethermalization plateau is reached. A conser¬ 
vative bound for the allowed degree of disorder can be es¬ 
timated from the presented linear response analysis, how¬ 
ever, a more realistic calculation must take into account 
self-regulation and slowing down of the unstable modes in 
the presence of disorder, which is a computationally chal¬ 
lenging task. On related grounds, it is desirable to study 
the formation of topological defects in quenches to the 
ordered phase, corresponding to the instantaneous limit 
of the quantum Kibble-Zurek mechanism [77, 78]. For 
the 3D Heisenberg model with SU{2) symmetry, topo¬ 
logically stable hedgehogs [79] are expected to form with 
universal scaling laws. Other possible research directions 
include extension to open spin systems, studying the role 
of NNLO corrections to assess the robustness of the NLO 
results, and comparison with other systematic expansions 
such as 1/D-expansion in D-dimensional lattices, and the 
semi-classical l/5'-expansion. 
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Appendix A: Summary of the truncated Spin-2PI 
formalism at the NLO level 


In this Appendix, we provide supplementary mate¬ 
rial for the Spin-2PI formalism along with a brief ac¬ 
count of the numerical methods. The covered material 
includes the explicit derivation of the approximate dy¬ 
namical equations from the NLO truncated 2PI effective 
action and the reconstruction of real-time spin-spin cor¬ 
relators from the Bethe-Salpeter equation. 


1. Correlation functions on the Schwinger-Keldysh 
time contour 


In the Schwinger-Keldysh formalism, the non¬ 
equilibrium dynamics of quantum fields is most elegantly 
derived from a path-integral defined on the round-trip 
contour C = U C“: 


C+ 

t = to \ t = +CXD 

c- 

The Majorana operators r] and real vector boson (p are 
replaced by Grassmann and real vector valued variables 
in the path-integral, along with an anti-periodic and pe¬ 
riodic boundary condition at the contour endpoints, re¬ 
spectively. 

The correlation functions defined on the C contour can 
be thought of 2 x 2 matrices in the two-dimensional space 
of the contour branch index. Eor example, the Majorana 
2-point correlator Q can be explicitly written as: 






(Al) 


where the times appearing in the matrix are ordinary 
times. We have dropped the discrete indices for brevity. 
The off-diagonal matrix elements are identified with the 
“lesser” and “greater” explicitly ordered correlators: 

G~~'~{ti,t 2 ) = G^{ti,t 2 ) = -i{r]{tx) r]{t 2 )). (A2) 

The diagonal matrix elements are related to each other 
by the virtue of the unitarity of evolution: 

G^^{tl,t2) = +0{ti - t2)[G^ {tl,t2) - G^{tl,t2)\, 

G (il,i 2 ) = -0{t2 - ti)[g> (ti,t2) -G^{ti,t2)\, 

(A3) 


which are identified with the usual retarded and advanced 
response functions, and Q While the 

lesser and greater correlation functions are independent 
functions for Dirac (complex) fermions, they are related 
to each other for Majorana fermions by transposition and 
negation, as it can be seen from Eq. (A3): 

g>(l,2) = -C?<(2,l). (A4) 


In summary, the 2-point correlator of Majorana fermions 
on the contour is fully specified by a single real-time 
correlator, e.g. ^>(1,2). It is easily shown that the 
same decomposition and relations hold for the Majorana 
self-energy U. The correlator of real bosons V and the 
bosonic self-energy U admit a similar decomposition, ex¬ 
cept for the absence of the relative minus sign in the 
definition of and : 


P+ {ti,t2) = T>‘^{ti,t2) = -i{(p{t2)(p{ti)), 
T^~~''{ti,t2) = T>>{ti,t2) = -i{(p{ti)(p{t2)), (A5) 






12 


which imply: 

P>(1,2) =P<(2,1). (A6) 

Similar to Majorana correlators, the 2-point correlator of 
real bosons on the contour is fully specified by a single 
real-time correlator, e.g. . The same result holds for 
the bosonic self-energy 71. 

2. Feynman rules for the Spin-2PI formalism 

The conventional Feynman diagram rules are used 
for interpreting the diagrams appearing throughout this 
work: 

i^ww>^2 zD“/5(1,2) 

l!-f2 ^ 

The integer indices refer to the bundle of lattice site and 
contour time in the diagrams above. Since the Majorana 
fermion propagators possess no charge flow direction, one 
may arbitrarily assign a direction to each line. The over¬ 
all sign of each diagram, however, must be determined at 
the end by counting the number of fermion permutations. 

The power counting of the large-extension is per¬ 
formed as follows: (1) each Majorana fermion loop intro¬ 
duces a factor of N resulting from the replica summation, 
(2) each interaction and boson line introduces a factor of 
l/N. 

The vacuum diagrams accompany symmetry factors 
which must be worked out case by case. The self-energy 
77 and 4-point vertex diagrams have an extra 
factor of i and respectively. 

3. Evolution of correlations functions in the 

Spin-2PI formalism 

The transition from the path-integral to the 2PI ef¬ 
fective action V[Q^V^Lp\ was briefly outlined in the main 
text and is a straightforward generalization of the results 
of Cornwall, Jackiw, and Tomboulis [41]. Within this for¬ 
malism, the evolution equations follow from a variational 
principle, reminiscent of Lagrangian dynamics of classi¬ 
cal particles, with the quantum correlators playing the 
role of generalized coordinates. Going back to Eqs. (8a) 
and (8b) and making V stationary with respect to P, 
and (^, we obtain: 

g-^ = g^^ - M[<pj - sig, v], (A7a) 

V-^ =Vo^ -n[g,v], (A7b) 

1 ^ 

i). (A7c) 


With the spin, replica, time, and space indices laid out 
explicitly, the “bare” fermion and boson propagators are 
written as: 

Oq (f 5 2) — ^aia2^aia2^jij2 ^ 

= N {V-X\T Sc{h,t 2 ), (A8) 

respectively. The contour Dirac (^-function is defined as 
— t 2 ) with the ± sign corresponding 
to ti,t 2 G C^, respectively. In Eq. (A7a), M[(^^](l,2) = 

the LO interaction 

effect and describes the coupling of Majorana fermions 
with the classical spin mean-field According to 

Eq. (A7c), the latter is instantaneously determined by 
the Majorana tadpole contracted with a bare interaction 
line. Thus, we find: 

^0-1(72 ^C(^l 5 ^J) ^^^7^ 

1 ^ 1| 

^ E = 2 X , (A9) 

cr = l 

which resembles the familiar Hartree self-energy that 
describes the mean-field effects. We emphasize that the 
Majorana tadpole is identified with the magnetization 
in our formalism. Note that the M(l,2) oc (^c(^i 5 ^J) is 
instantaneous and carries no memory effect. We will 
later show that truncating the approximation at this 
level and neglecting the self-energies indeed yields the 
Bloch equation. 

Going beyond the LO approximation, U[Q^V] and 
n[Q^V] describe memory effects associated from the spa- 
tiotemporal fluctuations of the exchange field. By defini¬ 
tion, these self-energies are obtained from the variations 
of r 2 [ 0 ,P]: 

/r|5,D|(l,2)s22i|^. (AlO) 

We recall that r2[0,P] is formally equivalent to the sum 
of 2PI vacuum diagrams constructed from the interaction 
vertex Cint[r],(p] = ^ ^a^'y admits a sys¬ 

tematic expansion in l/N. Here, we truncate the series 
at the NLO level: 

= ltr[PiTo] (All) 

where 77^(1,2) = «e^a/3 Ef=i y“y2’^’'^(772) 

gflX ’^is the Majorana bubble. Since V ~ l/N 
and the factor of N resulting from the replica summa¬ 
tion in the Majorana bubble, we find ^ C^(l)- This 

must be compared to the LO term in V int which is 0{N). 
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The resulting NLO self-energies are given as: 

2) = 4 X 1 

^NLO(i 2)^2 x iv^^^^^^^ = l77o(l,2). (A12) 

Having derived the explicit expressions for the self¬ 


energies, we discuss the derivation of evolution equa¬ 
tions as the next step. Our starting point are the cou¬ 
pled Dyson’s equations given in Eqs. (A7a) and (A7b). 
Strictly speaking, Dyson’s equations are differential iden¬ 
tities on the contour functions. They can be cast into a 
more useful form by acting them from the left and right 
hand side by Q and P, respectively, resulting in a set of 
contour integro-differential equations: 


^kj2^t2 ~ 


S{1,2) + J^dr rjy (i 1 , r) (r, ^ 2 ), 

<5(1, ^ ^hk (^1’ '^) * 2 ), 


(A13a) 

(A13b) 


(ti,t2) = ^y:r'^c(ti,t2) + ^y“r / dr77,y(ta,r)D-=(r,i2), (A14a) 

J c 

= + ^ / drD“y(ta,r)i7'fr(r,i2)V^r- (^Wb) 

J c 


We have defined shorthand ^(1, 2) = ^aia 2 ^ 2 ) 

and the contour integral dt A{t) is interpreted as 
G C+)-/^^dt^(t G C-). Eq. (A13a) 
and its adjoint Eq. (Al3b) are referred to as Kadanoff- 
Baym (KB) equations. The convolution integrals of self¬ 
energies and correlators manifestly show memory effects, 
which is a shared feature of beyond mean-field approxi¬ 
mations. 

The spatial structure of Eqs. (A13a)-(Al4b) can be 
simplified by noting that physical initial states imply ini¬ 
tial correlations between pairs of Majorana operators on 
the same site, i.e. (^oAo) oc Sj^j^. Crucially, this 

property extends to all times in the KB dynamics, inde¬ 
pendent of the order of truncation in 1/A^. To see this, 
one first establishes that the assumption Qj^j 2 {^iA 2 ) oc 
^iii 2 for to < ti,t 2 <T implies Uj^j^{tiA 2 ) oc 5j^j^ in 
the same domain. The causal structure of Eqs. (A13a)- 
(Al3b) subsequently extends this property to an in¬ 
finitesimally larger domains, and eventually to all times 
by induction. Therefore, we can always make the follow¬ 
ing simplifying substitution in the KB equation: 

^ (A15) 

The LO approximation: The KB equations reduce to the 
mean-field Bloch equation upon truncation at the LO 
level which amounts to neglecting fluctuation self-energy 
corrections —> 0. In this limit, the KB equations imply: 


-idt^Gfp'"’^{h,t2)+i(Pcj{t2)e^^,c,2 = 0 . 

(A16) 

Subtracting the equations from one another, setting t 2 = 
ti = t and using Eq. (7), we finally obtain: 

dt{Sj{t))=cp^jX{Sj{t)), (A17) 

which is the Bloch equation as anticipated. 


The NLO approximation: Including self-energy correc¬ 
tions, the time convolutions appearing in the KB equa¬ 
tions prohibit us from arriving at a closed equation for 
the equal-time Green’s functions and we inevitably need 
to solve for the complete unequal time Green’s function. 
Eor concreteness, we consider the case of spin spirals 
hereafter. The spatial structure of the KB equations 
can be significantly simplified by applying the unwinding 
unitary transformation, either directly on Eqs. (A13a)- 
(Al4b) or on the spin Hamiltonian. Either way, the 
initial spiral state transforms into an uncorrelated x- 
polarized EM state |l^o) = I ~^)j expense of 

an anisotropic interaction (see Eq. 13). The 2-point cor¬ 
relator of Majorana fermions at t = to is easily found 
as: 


'^aia2 ,> 
^3132 


{to,to) = ( 5 , 


/-il2 0 0 -i/2\ 

0 -i/2 1/2 0 

0 -1/2 -i/2 0 

V-i/2 0 0 -i/2/ 

(A18) 

The exchange field correlator at t = to is not an indepen¬ 
dent degree of freedom and is determined by 0(to,to), 
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see Eq. (A12). For translationally invariant initial states 
as such, Q and U further become independent of the lat¬ 
tice site. Furthermore, depends only on the dis¬ 
tance and at the NLO level, is local as well. The 

simplified structure of the correlators and self-energies is 
summarized as follows: 

0(1, 2) ^ <^RiR2 ^ 2 ), 

17(1,2) 


77(1,2) 7T“^(ii,i2). 


The KB equations can be written explicitly in terms 
of 0^ and using the Langreth rules [80]. We quote 
the final result, setting N = 1\ 


Jto 

- [ ' dri:-^^’>(ti,r) [0^"->(r,t2) + 0"^^’>(t2,r)] , 
Jto 

-iat.0"^"Hti,t2)+i0"^^’>(ti,t2)(/^rH^2)= / 'dr [0-^^’>(ti,r) + 0^-^’>(r,t)] r^-^’>(r,t2) 

Jto 

- r dr0^^^’>(ti,T) [i:^^^’>(r,t2) + i:^^^’>(t2,T)] 
Jto 


^OiiCX2 


^aia2 


/ 'dr [71^"’>(ti,r) -71"^’>(r,ti)] {tM) 

Jto 

- r dT7T'‘"’>(7i,r) [pr'’>(r,i2) -7>r"’^(i2,r)] , 

Jto 

= F“^'‘7T^"’>(7i, 72) + / Tr [D“^'‘’>(7 i,t) - D'(“^’>(r,7i)] 77'‘"->(r,72) 

Jto 

- [ ' dTV^^^’>{h,T) [7T^'^’>(r,72) - 77‘^^’>(i2,r)] 

Jto 


(A19a) 


(A19b) 


(A20a) 


(A20b) 


The self-energies and 77^ are read from Eq. (A12). 
The last explicit equations are suitable for devising a 
numerical forward propagation scheme. Starting from 
00 (^ 05 ^ 0 ): we calculate i7(to,to) and 71 (to, to) from 
Eq. (A12), and P(to,to) from Eq. (A20a). The casual 
structure Eqs. (A19a)-(A20b) allows us to propagate 
{0,i7,P, 71} in (ti,t 2 ) in discrete time steps of size 
At. This is achieved using a robust predictor-corrector 
method with guaranteed accuracy to 0{At^). 

Calculating spin-spin correlators in the Spin-2PI formal¬ 
ism: In the framework of 2PI effective actions, higher 
order correlators are “reconstructed” from the history of 
2-point correlators. Here, we are interested in the con¬ 
nected spin-spin correlator: 

= (n[sy(<i)s«fe)]) 

-(»?(*■)) (s"fe)). (A21) 

where the spin operators are shorthand notations for 


Eq. (4). The spin-spin correlator is found from the Ma- 
jorana L-function, defined as: 

L(ll; 22) = {Tc [v{l) v{l) g{2) g{2)]) - ig{C 1) ^0(2, 2), 

(A22) 

by contracting 5-symbols with its left and right pair of 
fermion lines: 


X 


■3132 


'{tiM) 


i 

^ ^ai/3i7i 


r/3i7i;/3272 

3131-^3232 


(t^ , ti ; , ^2 ) ^^2/3272 • 


(A23) 

The L-function in turn satisfies a non-equilibrium Bethe- 
Salpeter equation on the C contour: 


7i(ll; 22) = 712(11; 22) 

+ J d3 d3 d4d4772(11; 33) (33; 44) 77(44; 22), 

^ (A24) 

where 772(11; 22) = 0(12)0(21) - 0_(12)0(l2) 

and the 2PI irreducible vertex (33; 44) = 

(5^rint[^7]/(50(33)(50(44); here, rint[0] is given in Eq. (8b) 
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with and V substituted in terms of Q from the 
stationarity condition Eqs. (A7b) and (A7c). 

The NLO effective action yields three contributions to 
A( 2 ): 


^( 2 ) 




AAAAA/T—4 


+ 8 x 


2_iAAAAA/i_4 


A 


( 2 ) 

MT 


A 


( 2 ) 

AL 


+ 



(A25) 

The last symbol stands for the three permutations of 
the first three diagrams obtained by (3 3), (4 4), 

and (3 3,4 4) with signs —, — and +, respec¬ 

tively. These vertex corrections are structurally simi¬ 
lar to the RPA, Maki-Thompson (MT), and Aslamazov- 
Larkin (AL) vertex corrections accounting for supercon¬ 
ducting fluctuations in metals [81]. Explicitly, these ver¬ 
tex parts are given as: 


^rpa(^^5 44) — 


44) — 


X Sc{t3:t^)Sc{U,ti) Scits^U), (A26a 

X 5c{h,U)5e{h,ti), 


*^5175(^35^3) r, Sj 3 j 4 


(A26b) 


4'^(33; 44) = {ts, ts) (^ 4 , h) 


J 3 J 3 

X 2 (^3: ^4) “ 

X - (^3: ^i) A ^vcx^l3^ 


(A26c) 


It is easily noticed that ~ ^ 0(1/A^), 

and A^l ^ 0(1/A^^). Therefore, we may drop the latter 
if accuracy at the NLO order is desired. 

Eor translation invariant states, it can be shown that 
L(n;22) - Taking 

a Eourier transform in Ri — R 2 yields decoupled inte¬ 
gral equations for each momentum transfer q. The tem¬ 
poral structure of the BSE remains formidable. Per¬ 
forming the contour integral and discrete summations 
over 3,3 variables in Eq. (A24) and contracting the 
right legs according to Eq. (A23), we reach to an in¬ 
tegral equation for the 3-time object Tqtj;^ 2 ) = 
{-i/ 2 )Lq^°‘^'°‘^°‘^{ti,ti;t 2 ,t 2 )Sf^a 2 a 2 in two contour 
times (with fixed external time ^ 2 ). 

The first step in solving the BS equation is to recast 
it in terms of functions of ordinary times. In compari¬ 
son to the KB equation, this step is significantly more 
involved here due to the complex real-time structure of 
3-time and 4-time CTP functions and multiple contour 
integrals. We leave the cumbersome details for a sepa¬ 
rate publication and solely outline the procedure here. 
We showed earlier in Sec. A 1 that the 4 real-time matrix 
elements of 2-time functions such as Q and V can be fully 
specified using a single real-time function, e.g. . A 
similar analysis of Tq tj; ^ 2 ), taking into account 

symmetries and unitarity of evolution, reveals that the 8 
real-time components of 3-time function as such can be 
fully specified by 3 independent functions. Accordingly, 
the contour BS equation can be explicitly written as 3 
coupled two-dimensional integral equations in ordinary 
times; the latter is numerically solved by discretizing the 
integrals using approximate quadratures and solving the 
resulting linear system. 
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